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Abstract 

t I Solid-On-Solid (SOS) computer simulations are employed to investigate 

the sublimation of surfaces. We distinguish three sublimation regimes: layer- 
by-layer sublimation, free step flow and hindered step flow. The sublimation 
regime is selected by the morphology i.e. the terrace width. To each regime 
' corresponds another effective energy. We propose a systematic way to derive 

microscopic parameters from effective energies and apply this microscopical 
analysis to the layer-by-layer and the free step flow regime. We adopt analyt- 
ical calculations from Pimpinelli and Villain and apply them to our model. 
■ keywords: Computer simulations; Models of surface kinetics; Evaporation 

and Sublimation; Growth; Surface Diffusion; Surface structure, morphology, 
roughness, and topography; Cadmium telluride 
PACS: 81.10.Aj; 68.35.Fx; 68.10.Jy 



1 Introduction 



The sublimation of semiconductor materials such as cadmium telluride or silicium 
has been of considerable interest in the last few years. The motivation of these 
studies has been mainly to understand the kinetics on these surfaces. This would 
be helpful e.g. to gain better control in the fabrication of nanostructure devices. In 
this paper we want to describe computer simulations of sublimation. We use the 
conceptually simple SOS model which has been extensively used to model Molecular 
Beam Epitaxy (MBE) [|I|, E3. As far as we know these are the first simulations of 
sublimation of this kind. 

Theoretical studies of sublimation were based on the 1+1 dimensional Burton- 
Cabrera- Frank model. In this model the evolution of the surface is described by the 
motion of steps. In 2+1 dimensions it can be interpreted as a surface of parallel, 
mono-atomic and straight steps. Pagonabarraga et. al. examined the validity of 
the continuum approximation in the step flow regime |§. Pimpinelli and Villain 
discussed the creation of surface vacancies and found a condition for the occurrence 
of "Lochkeime" (poly- vacancies) |4[] . In the following we will refer to calculations in 
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their introductory book of crystal growth ||. In particular they found that in step 
flow sublimation the velocity of steps is influenced by the terrace width. 

Experimental studies of CdTe(OOl) concentrated on the extraction of effective en- 
ergies. Several authors used Reflection High Energy Electron Diffraction (RHEED) 
to study sublimation || [7|, ||. Other investigations were based on the use of a 
quadrupole mass spectrometer (QMS) |5|, [TU| (and references therein). Recently 
Neureiter et al applied Spot Profile Analysis of Low Energy Electron Diffraction 
(SPA-LEED) to CdTe(OOl) [0, [|. The oscillation period r in RHEED or SPA-LEED 
intensity can be described by an arrhenius-law 1/r = v e~ E l kT '. For CdTe(OOl) an 
activation energy of l.QeV - derived from the diffraction data - is commonly accepted 
. Using QMS a value of 1.55eV^ for the mass desorption rate has been reported. We 
will show that this discrepancy can be explained by the existence of different mor- 
phology dependent sublimation regimes. Neureiter et. al. show that such a picture 
of the CdTe(OOl) surface is consistent with SPA-LEED and QMS measurement s|9|]. 

It would be desirable to obtain microscopic parameters of a computer model 
from ab-initio calculations. Calculations are now available for GaAs(OOl) [12|, [13 



but for II- VI semiconductors no such calculations are available yet. Thous computer 
simulations of MBE or sublimation can help to get insight into these parameters. 
E.g. Smilauer and Vvedensky were able to find a set of parameters for GaAs(OOl). 
They compared the step density of their model to RHEED intensities during MBE 
growth and subsequent recovery |^4| . These investigations showed that it is possible 
to describe such a complicated surface (reconstruction, 2 types of atoms ...) by 
means of the simplest SOS model. Naturally the microscopic parameters found are 
still effective parameters. We expect that these values should be related to the true 
microscopic interactions describing the slowest processes. 

There is a certain advantage in using sublimation data compared to MBE to 
find microscopic parameters of real materials. First of all the process is determined 
uniquely by microscopic processes, whereas in MBE the evolution of the morphology 
is dominated by the external flux of atoms. Secondly, sublimation has been recently 



studied by SPA-LEED [11, 0[. To a great extent this method can be analyzed in 



kinematic approximation [JT3L [111, whereas RHEED-profiles have to be calculated 



using multiple-scattering theory [pL7| . Also the mass desorption during sublimation 
can be easily compared with computer simulations. 

The organization of this paper is as follows. In section two we will introduce 
the SOS model and the algorithm used in the simulations. In section three we 
will investigate sublimation similar to experimental investigations. We will identify 
three sublimation regimes in dependence of the terrace width. In section four we will 
connect effective activation energies (which are measurable in real experiments) to 
the underlying microscopic parameters. On the one hand we will be able to compare 
these relations with theoretical predictions of Pimpinelli and Villain (section five). 
On the other hand this leads to a systematic approach for determining microscopic 
energies of real materials which will be illustrated in section six. A discussion and 
an outlook will follow in section seven. 
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2 Solid-On-Solid model 



One basic assumption in this model is the neglection of overhangs. Thus no bulk 
vacancy formation can be investigated. The second assumption is a fixed lattice 
structure. Both lead to the crucial simplification for the simulation that the surface 
can be stored in one simple data-array of heights. 



Figure 1: The activation 
energy in the jump rate of 
atoms depends on the local 
arrangement. The diffusion 
barrier is Eb, for desorption 
an energy of E D is neces- 
sary. Next in-plane neigh- 
bours have a binding en- 
ergy of E N . A jump over 
a step edge is hindered by 
the Ehrlich-Schwoebel bar- 
rier E s . 



We have used a simple version of this model: a simple cubic lattice with only one 
species of atoms. Since there is only one lattice constant all lengths will be given 
in units of this lattice constant. Figure [I] shows some selected processes and their 
corresponding activation energies. 

Theories of activated processes lead to an arrhenius form of the jump frequency 

-E/kT 




z/ e 



The prefactor in general depends on the form of the potential and 
the temperature. Nevertheless the exponential factor gives the dominant tempera- 
ture dependence. Thus we assume a constant prefactor z/ . The activation energy 
depends on the microscopic mechanism involved. We assume that diffusion has an 
activation energy of E = Eb + n ■ E^ where n is the number of next in-plane neigh- 
bours. Thus free atoms in this model have the diffusion constant D = uq e~ Es ^ kT . 
At step edges an Ehrlich-Schwoebel barrier Eg is considered. In addition to "con- 
ventional" SOS-models of MBE we take into account the desorption similar to the 
diffusion with an jump frequency T desorption = u e ~( E D+n-E N )/kT _ ^y e assume fo e 
same prefactor for desorption as for diffusion which is certainly not valid in real 
materials. 

Especially when simulating sublimation it is absolutely important to include all 
possible events because even events with very low rates can become very important. 
E.g. the creation of a surface vacancy has a rate which is several orders of magnitudes 
lower then the jump frequency of free adatoms. Nevertheless this process is the most 
important one to initiate layer-by-layer sublimation. 

The Maksym-algorithm JL9J is able to handle such big differences in jump rates. 
In this algorithm all possible events are considered simultaneously and are stored in 
one array of rates. The event with a high rate will be chosen with a corresponding 
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high probability. Maksym described how to find in an efficient way the event selected 
by a random number in the range (0, grates). We improved the simulation of the 
(time consuming) diffusion of free adatoms by handling them separately. 

In this paper we will concentrate on two sets of parameters: Vq = 10 s , Eb = 
0.9eV, E D = l.leV, E N = 0.25eV, and E s = O.leV (set 1). A second set u = 
10 12 s-\ E B = 0.5eV, E D = 0.7eV, E N = 0.35eV, and E s = O.leV (set 2) will be 
compared. 

3 Crossover from layer- by- layer to step flow sub- 
limation 

Since computer simulations are restricted to rather small systems (we use system 
sizes up to 800 x 800 lattice constants, typically 512 x 512), one is forced to induce 
equidistant steps by fixing a height difference between two opposite boundaries. We 
have investigated the effect of such an induced step train on the mass desorption. 

There are two limiting cases. First the pure layer-by-layer sublimation. In this 
case the poly-vacancies are much smaller then typical terrace sizes. Changing the 
terrace sizes in this regime should not influence the total amount of layers desorbed 
after a fixed time. Furthermore the number of layers desorbed should be equal to 
the number of periods observed in LEED or RHEED. 

In the second extreme the steps are so close together that there is no place to 
create poly-vacancies and there will be only desorption from steps. In this case the 
number of desorbed monolayers Ah after the time At will be equal to 

Ah _ V s tep 

At Lf errace 

because L terrace /v step is just the time to evaporate one monolayer. 



Figure 2: Height loss per 
time over 1/L terrace . Each 
data point represents a sim- 
ulation of parameter set (1) 
at 600K. For all simulations 
At was chosen to be 800s. 
The dashed line corresponds 
to v step = 0.154s" 1 . 
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Figure |2| clearly show the existence of these two regimes. However at very small 
terrace sizes a third regime can be identified (we have not plotted another data 
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point at 1/ Lterrace = 0.2 which clearly establishes region III). This behaviour was 
already calculated by Pimpinelli and Villain using the modified Burton- Cabrera- 
Frank (BCF) model||. In this framework but without a Schwoebel-effect (however 
the estimation is valid in general) v step oc tanh(L termce /2 y/Dro) , where D is the 
diffusion constant of free atoms and td the life-time of an particle until desorption 
occurs. If the terrace width is greater then the diffusion length ^/Dtd one gets a 
constant velocity. But if one lowers L terrace the particles will cross the terraces and 
will be reincorporated at the step edges. In this case v step will be proportional to 
Lterrace and Ah/ At will become constant again. That is why we want to distinguish 
between hindered and free step flow. Setting the argument of the tanh to one and 
identifying D = uq e~ EB ^ kT and 1/td = vq e~ ED ' kT one can estimate the crossover 
length between regime II and III to be £ x ~ 14 in good correspondence to figure 

Likewise Pimpinelli and Villain calculated the crossover between the free step 
flow and the layer-by- layer regime. In the approximation described in section five 
we obtain for parameter set (1) a crossover length of about 60, which is once again 
in good correspondence to our simulations. However we will show in section five 
that this agreement of theory and simulation is merely by chance. 

In addition we have investigated the intensity in antiphase of the (0,0)-spot using 
the kinematic approximation. We found that the oscillation period was the same 
either by calculating the step density or the width of the height distribution (in the 
case of singular surfaces). The evaluation of the LEED-profiles is complicated by 
the fact that due to the small system sizes the distribution of terrace sizes is limited. 
This leads to a pronounced twin-peak structure in k-space at k vea k = i-p-^ 1 — . We 
have mimicked the effect of a broad distribution of terrace sizes by convolution of the 
obtained profile with an Gaussian of 0.5% of the Brioullin Zone. In correspondence 
to the conclusions drawn from the mass desorption diagram (fig. Q) one oscillation 
period of 232s±8 down to a terrace size of Lterrace ~ 100 is observed. At Lterrace = 64 
an oscillation could be identified in the twin peak structure. However the period 
could only be estimated to be roughly r « 265s ± 20s. With L terrace = 32 no such 
oscillations occurred. 

We have shown that in regime I (layer-by-layer sublimation) the observed oscil- 
lations are independent of the terrace size. Decreasing the mean terrace width one 
reaches the free step flow regime (regime II) which is characterized by a constant 
step velocity. Furthermore oscillations do not exist in this regime. However at small 
terrace sizes (small compared to the diffusion length) one reaches the hindered step 
flow regime III which is again characterized by a mass desorption rate independent 
of the terrace width. 

In each regime one could expect to define effective energies. In regime I the 
oscillation period is a well defined variable. In regime II -the free step flow regime- 
the step velocity and at small terrace sizes (regime III) the mass desorption rate 
can be analyzed. We concentrate on regime I and II but we will mention the 
theoretical result (including the Ehrlich-Schwoebel-barrier) in section 5.3 for regime 
III. Looking at the temperature dependence of 1/r and v step (figure ||) one can 
indeed identify arrhenius laws with different activation energies. 
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Figure 3: Comparison 
of the activation energies 
of the layer-by-layer and 
the free step flow sublima- 
tion of parameter set (1). 
In the layer- by- layer regime 
[Lterrace = 256) the oscilla- 
tion period is characterized 
by the effective activation 
energy Ej = 1.73eV. At 
Lterrace = 16 the step ve- 
locity yields E n = 1.60eV. 

4 Interpretation of effective energies: microscop- 
ical analysis 

We will now try to establish a link between the microscopic parameters of the model 
to the effective energies. The effective energies are derived from the variation of 
the temperature. In the simulations however it is possible to vary the microscopic 
parameters as well. 

We will vary the energies of set (1) by fixing three energies and changing the 
fourth. In detail we have set the parameters to E B /eV={0.8, 0.85, 0.9, 0.95}, 
E D /eV={1.00, 1.05, 1.10, 1.15}, E N /eV={0.23, 0.24, 0.25, 0.26} and E s /eV={0.05, 
0.10, 0.15, 0.20}. Assuming a linear dependence of the effective energy on the 
parameters -^regime = P^b + 8Ed + vE^ + &Eg we want to extract the importance 
of the different microscopic parameters. 

4.1 Layer- by-layer sublimation (regime I) 

Our simulations at different temperatures (s. figure |3| ) show an arrhenius behaviour 
with vf = 1.48 10 12 s _1 and an effective activation energy of Ej = 1.73 eV. The 
upper index "T" indicates that these are values obtained by simulations at different 
temperatures. 

The investigation of the effect of the microscopic parameters yields v l {°" = 6.3 10 11 s 
and 

E l j a - = 0.61 ■ E B + 0.35 • E D + 2.85 ■ E N + 0.44 • E s . 

In our case this leads to an effective activation energy close to the "conventional" 
effective energy 

E l j a - = l.Q9eV ^Ej = 1.73eV. 

Note that E\°" has been deduced from simulations at a fixed temperature T = 600^. 
At first glance it is surprising that the diffusion parameter Eb has a larger influence 
on the effective energy than the desorption energy Ed- This can be explained by 
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the typical life of an evaporating particle. Firstly the atoms become freely diffusing 
atoms (this needs an activation energy of typically Eg + 3- E^. Thereafter the atoms 
diffuse on the terrace until they evaporate. The Ehrlich-Schwoebel barrier plays an 
essential role in the creation of surface vacancies which explains the influence of Es- 
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Figure 4: Temperature de- 
pendent sublimation rate in 
the layer- by- layer regime of 
set (2). However the pre- 
dicted curve was multiplied 
by a factor 1.2 to show that 
the predicted energy is very 
good. 
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Table 1: Predicted and measured Oscillation periods. The simulations were carried 
out at T = 600i^, the vibration frequency was set to z/ = 12 12 s _1 

Other simulations show, that this approximation of first order is valid over a 
greater range in parameter space. The predicted periods of these models are com- 
pared with measured ones in Table |I| Figure f| shows that the activation energy 
predicted by our linear approximation is still very good for the set (2) although the 
prefactor no longer holds. 

4.2 Free step flow (regime II) 

In the following investigations we used a surface with a mean terrace size of 16 
lattice constants. Looking at Figure |2| one can see that one is in the free step flow 
regime. 

By evaluating the temperature dependence in this extreme limit we find ujj = 
3.9 • 10 12 s- 1 and E T U = l.QOeV. 

As in the case of layer-by-layer sublimation we connected this effective relation 
to the microscopic parameters. Of course we cannot be sure to stay in the free step 
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flow regime. We checked that Ah/ At is at least twice as high as in the pure layer- 
by-layer regime. However with high Ed or low Eb we reached the hindered step 
flow regime due to the fact that the atoms stayed to short a time on the terrace. 
Our results can be fitted by 



E 



La. 
U 



0.70 • E B + 0.32 • E D + 2.4 • E N + 0.21 • E s 



1.60eV. 



In this case the two different effective energies are identical Ej£ = 
Likewise the prefactor of 4.3 • 10 12 s^ 1 is close to the previous one. 

Comparing this relation with the equivalent relation for the layer-by-layer regime 
one observe that the contribution of E^ and E$ have significantly changed. This 
seems plausible because the creation of a surface vacancy compared to the creation 
of a kink is hindered by the energy En ~\~ Eg. 
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Figure 5: Predictions and 
measurements of v step for 
parameter set (2). No cor- 
rections were necessary! 
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A comparison between our prediction for the set (2) and measurements are shown 
in figure |||. As one can see our linear approximation is very good in this case. 

5 Comparison with theoretical calculations 

Now we are able to compare the theoretical predictions of Villain and Pimpinelli 
with our results. We will interpret variables in their results in terms of arrhenius 
activated jump rates. 



5.1 Free step flow (regime II) 

In the free step flow regime they obtained the following expression for v s t ep on a 
surface with terraces of equal size £ (ch. 6.4 in 0) 



J step 



PqkD 



4 - (2 - kD (±r + i)) e-** - (2 + kD + ^)) e 
(l + k§) (l + K§r) 



1 - K 



1 — K ~§r^j e 



8 



Po can be identified with the ratio between the detachment rate from a step and 
the diffusion jump rate (p = e~( n ' EN ^ kT , where n ■ E N is a typical binding-energy 
to next in-plane neighbours at a step). The diffusion parameters at step edges are 
related with the diffusion constant D. Diffusion downward a step D' = e~ Es ^ kT D. 
Diffusion towards a step from the lower terrace is not changed in our model (D" = 
D = i/ • e~ EB / kT ) . The inverse of the diffusion length is given by k — 1/ \JDtd = 
e -(E D -E B )/2kT _ Ta k ing t j ie limit for nl > 1 (free step flow) and neglecting e -( E D-E B )/kT 
versus 1 we obtain the simpler relation 

E B /2+E D /2+n-E N / 1 \ 

^ tep ~ ^ • ^1 + x - e _ {ED , 2 _ EB , 2 _ Es) ,Kr ) ■ 

This shows that the assumption of an arrhenius behaviour is not valid in this 
sublimation regime. However with parameter set (1) the correction factor has 
no temperature dependence . Because E D /2 — E B /2 = E s we obtain v step ~ 
3/2 i/ Q Q-^ E B/ 2 + E D/^+n-E N )/kT _ j n microscopical analysis we measure an effec- 
tive importance of E D , E B and E s due the correction factor. The influence of 
E N is unchanged. Indeed if we use the theoretical result with n ■ E n = 2.4 • E N 
(as the microscopical analysis suggests) we obtain an effective energy of 1.6eV in 
correspondence to our measurement. 



5.2 Crossover to layer- by- layer sublimation 

Given the two equations for the period of oscillations and the step velocity we 
can calculate in our linear approximation approach the dependence of L x from the 
temperature and the microscopic parameters. From 1/r = v ste p/L x we obtain 

L _ u fsf c -(E f „f-E m )/kT ^ g g e (0.03£ D -0.09£ s +0.45£jv+0.23£ s )/fcT 

m 

We already mentioned the theoretical relation for the crossover terrace width. Pimpinelli 
and Villain estimated the crossover to be at 

\ 2 n /\ , Aff o\ 



,2 x ln£ x ^4tt (jt-) t £>D 1 1 



kTJ \ DpoJ 

(To is the equilibrium concentration and A the diffusion constant of the surface va- 
cancies. 7 is the free energy per step length. Neglecting A and setting 7 = E N 
(which is a good approximation if the step is nearly free of kinks) we can relate 
this expression to our model. If we simplify their relation in the vicinity of £ x at 
T = 600K we obtain approximately 

£x * ^ (§0 ■ e(ED ~ EB)/kT - 

which is at the same order of magnitude as our result but has a very different temper- 
ature dependence. In the theoretical result only desorption and diffusion energies 
contribute to the exponential temperature dependence. Whereas our simulations 
indicate that these terms cancel out and only E N and Es contribute exponentially. 
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5.3 Hindered step flow (regime III) 

We have not examined in detail the hindered step flow regime by computer simu- 
lations. Nevertheless we want to mention the theoretical result for this regime. In 
this regime rf < 1 and therefor one can use the Taylor series of the exponential 
functions e m 1 ± k£. In this case the step velocity is described by an arrhenius 
behaviour. 

Vste P « £u e~^ +n - E ^ kT 

This means that the number of monolayers desorbed is independent of the terrace 
width in this regime. 

6 Systematic derivation of microscopic energies 

Our connection between effective and microscopic energies will help to develop a set 
of parameters describing e.g. the CdTe(OOl) surface. Taking the effective energies 
for CdTe(OOl) as Ej = 1.9eV and Ejj = 1.55eVQ we can then reduce the four 
dimensional energy-parameter space to a two dimensional one. E. g. the Schwoebel- 
barrier can be expressed in terms of the diffusion and desorption barriers: Eg = 
0.31eV + 1.16E B + 0.16E D . In figure |6| the E N (E B , E D )-p\a.iae is shown. But because 
there are further physical restrictions (namely that E B > 0, E N > and E D > E B ) 
only a triangle is left from the infinite plane. One should however be careful. In 
the dark region (where Ed < E B + Eg) surface vacancies will be created via direct 
desorption rather than the creation of an adatom-vacancy pair. In this region our 
microscopic analysis is certainly incorrect. It should be emphasized that our relation 
of effective to microscopic energies is only valid in the vicinity of parameter set (1). 
However one can take it as a first approximation. Then one can systematically 
improve the result by iterating the microscopic analysis. 




Figure 6: Proposed param- 
eter plane E N (E B , E D ) for 
CdTe(OOl). Due to physical 
constraints a triangle is cut 
out. In the dark region our 
microscopical analysis is not 
valid. There the process of 
the generation of surface va- 
cancies is changed. 



1 It is possible that the mass desorption energy is due to regime III. We just want to demonstrate 
how the microscopical analysis can be used. 
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7 Conclusion 



In this paper we have reported the first simulations of sublimation in the framework 
of the SOS model. We have been able to identify three regimes of sublimation by 
looking at the mass desorption: layer-by- layer sublimation, free and hindered step 
flow. The existence of these regimes is in correspondence to predictions for the 
BCF-model from Pimpinelli and Villain. 

We have related the corresponding effective activation energies to the underlying 
microscopic parameters. Although our relation is just a linear approximation in 
the vicinity of the parameter set (1) the application to other sets of parameters is 
quite successful. From these relations we were able to compare the step velocity and 
the cross-over to layer-by-layer with analytical results from Pimpinelli and Villain. 
Our result for the step velocity is similar to the analytical result. However our 
approximation for the crossover length is very different from the analytical one. 

Our simulations show that the effective desorption energy is to a great extent 
influenced by the diffusion of adatoms. Thus energies measured in real experiments 
cannot be understood as desorption energies. Furthermore the energies measured 
are highly influenced by the morphology of the surface. We have only compared 
singular and stepped surfaces but it is now obvious that e.g. in MBE one has to 
expect other desorption rates than measured by Sublimation. 

In this work we have concentrated on the sublimation period and mass desorp- 
tion. Morphological features as typical surface vacancy distances will be the work 
of further investigations. This way we hope to get a more detailed insight in the 
physics of sublimating surfaces. 
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